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Abstract 

We discuss a novel analysis method for reaction network systems with poly¬ 
nomial or rational rate functions. This method is based on computing tropical 
equilibrations defined by the equality of at least two dominant monomials of op¬ 
posite signs in the differential equations of each dynamic variable. In algebraic 
geometry, the tropical equilibration problem is tantamount to finding tropical 
prevarieties, that are finite intersections of tropical hypersurfaces. Tropical equi¬ 
librations with the same set of dominant monomials define a branch or equiva¬ 
lence class. Minimal branches are particularly interesting as they describe the 
simplest states of the reaction network. We provide a method to compute the 
number of minimal branches and to find representative tropical equilibrations 
for each branch. 


1 Introduction 

Networks of chemical reactions are widely used in chemistry for modeling cataly¬ 
sis, combustion, chemical reactors, or in biology as models of signaling, metabolism, 
and gene regulation. Several mathematical methods were developed for anal¬ 
ysis of these models such as the study of multiplicity of steady state solutions 
and detection of bifurcations by stoichiometry analysis, deficiency theorems, 
reversibility, permanency, etc. [1]. 

All these methods focus on the number and the stability of the steady states 
of chemical networks. Beyond steady states, metastable states defined as regions 
of the phase space where the system has slow dynamics are also important for 
understanding the behavior of networks. For instance, low dimensional iner¬ 
tial or invariant manifolds gather slow degrees of freedom of the system and 
are important for model reduction. Invariant manifolds can lose local stabil¬ 
ity, which allow the trajectories to perform large, rapid phase space excursions 
before slowing down in a different place on the same or on another invariant 
manifold. [H]. Biological networks have often been modelled as discrete dy¬ 
namical systems, whose trajectories are sequences of states in a discrete space, 
see, for instance, the Boolean automata of Rene Thomas m- We think that 
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metastable states of continuous models are good candidates for representing 
states in a discrete model. It is, therefore, important to know how metastable 
states are dynamically connected. 

We showed elsewhere that tropical geometry methods can be used to approx¬ 
imate such invariant manifolds for systems of polynomial differential equations 
[nmnn]. The connection between original dynamics of a specific system 
pn] to such a method based on invariant manifold is described in [T3]. In a 
nutshell, the slowness of the dynamics on the invariant manifolds follows from 
the compensation of dominant forces acting on the system, represented as dom¬ 
inant monomials in the differential equations. We have called the equality of 
dominant monomials tropical equilibration [niin]. Tropical equilibrations are 
different from steady states, because in tropically equilibrated systems one has 
non-compensated weak forces that drive the system slowly, whereas at steady 
state, the net forces vanish. Furthermore, invariant manifolds can be roughly 
associated to metastable states because they are regions of phase space where 
systems dynamics is relatively slower. In this paper, we introduce methods to 
compute tropical equilibrations and group them into branches that cover the 
metastable states of the system. These branches of tropical equilibrations form 
a polyhedral complex. The zeroth homology group (or in other words, the num¬ 
ber of connected components) of this complex indicates the possible transitions 
between the metastable states. Additionally, we explore the structure between 
the equilibration solutions using directed graphs (to present the inclusion rela¬ 
tions among them) and undirected graphs (to present the connectivity among 
them). Furthermore, one of the biochemical applications of these equilibrations 
is constructing invariant manifolds and to that extent defining slow variables. 
We provide a way to identify such slow species and visualise them through 
heatmaps. Lastly, we benchmark our method against models obtained from the 
Biomodels database m and discuss the findings. 

2 Definitions and Settings 

We consider biochemical networks described by mass action kinetics 

^ ^ , l<i<n, (1) 
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where kj > 0 are kinetic constants, Sij are the entries of the stoichiometric 
matrix (uniformly bounded integers, \Sij\ < s, s is small), cy.j = (aj,..., a;!) 

are multi-indices, and .. .Xn". We consider that a{ are non-negative 

integers. At this point, we like to mention, there exist approaches to describe 
such a polynomial system in a graph theoretic way, i.e., by a weighted directed 
graph and a weighted bipartite graph and to study the number of positive solu¬ 
tions depending on the graph structure [S]. This approach uses decompositions 
of Newton polytopes to find that parts of the directed graph are related to the 
existence of positive steady state solutions. In our paper, we do not use these 
graph theoretic considerations and investigate the different problem of tropical 
equilibrations. 

In the case of slow/fast systems with polynomial dynamics such as Q, the 
slow invariant manifold is approximated by a system of polynomial equations for 


2 


the fast species. This crucial point allows us to find a connection with tropical 
geometry. We introduce now the terminology of tropical geometry needed for 
the presentation of our results, and refer to |12j for a comprehensive introduction 
to this field. 

Let /i, /2, ■ • ■, //c be multivariate polynomials, fi G C[a;i, a;2,..., a:„], repre¬ 
senting all or a part of the polynomials in the right hand side of 0. 

Let us now assume that the variables Xi, i G [l,n] are written as powers 
of a small positive parameter e, namely Xi = x^s^’, where Xi has order zero 
(there are Ci,di not depending on e such that 0 < Ci < Xt < di). The orders 
Qi indicate the order of magnitude of Xi. Because e was chosen small, are 
lower for larger absolute values of Xi. Furthermore, the order of magnitude of 
monomials is given by the dot product /i = {a, a), where a — (ai,..., a„). 
Again, smaller values of fi correspond to monomials with larger absolute values. 
For each multivariate polynomial /, we define the tropical hypersurface T{ f) as 
the set of vectors a G K" such that the minimum of {a, a) over all monomials 
in / is attained for at least two monomials in /. In other words, / has at least 
two dominating monomials. 

A tropical prevariety is defined as the intersection of a finite number of 
tropical hypersurfaces, namely r(/i, /2 ,..., fk) = f^i(^[i,k]T{fi). 

A tropical variety is the intersection of all tropical hypersurfaces in the ideal 
I generated by the polynomials /i, /2, ■ ■ ■ ,fk, T{I) = n/g/T(/). The tropical 
variety is within the tropical prevariety, but the reciprocal property is not always 
true. There exist algorithms to compute such tropical varieties as in m- 

For our purpose, we slightly modify the classical notion of tropical prevariety. 
A tropical equilibration is defined as a vector a G K" such that {a, a) attains 
its minimum at least twice for monomials of different signs, for each polynomial 
in the system /i,/2, • ■ ■,/fc- Thus, tropical equilibrations are subsets of the 
tropical prevariety. Our sign condition is needed because we are interested in 
approximating real strictly positive solutions of polynomial systems (the sum 
of several dominant monomials of the same sign have no real strictly positive 
roots). 

3 Branches of Tropical Equilibrations 

For chemical reaction networks with multiple timescales, it is reasonable to 
consider that kinetic parameters have different orders of magnitudes. 

We, therefore, assume that parameters of the kinetic models Q can be 
written as 

kj = kje'^G (2) 

The exponents 7^ are considered to be integer. For instance, the following 
approximation produces integer exponents: 

= round(log(fcj)/ log(e)), ( 3 ) 

where round stands for the closest integer (with half-integers rounded to even 
numbers). Without rounding to the closest integer, changing the parameter 
e should not introduce variations in the output of our method. Indeed, the 
tropical prevariety is independent of the choice of e. 

Of course, kinetic parameters are fixed. In contrast, the orders of the species 
vary in the concentration space and have to be calculated as solutions to the 
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tropical equilibration problem. To this aim, the network dynamics is first de¬ 
scribed by a rescaled ODE system 

(4) 

j 

where 

/ij(a) = - 1 - (a,aj), (5) 

and (,) stands for the dot product. 

The r.h.s. of each equation in 

Q is a sum of multivariate monomials in the concentrations. The orders fij 
indicate how large are these monomials, in absolute value. A monomial of order 
fj,j dominates another monomial of order ^ji if fj,j < 

The tropical equilibration problem consists in finding a vector a such that 

.min (7j + (a, otj)) = min (7^ -t (a, a.j)) (6) 

This system can be represented as a set of linear inequalities resulting into a 
convex polytope. The solutions of this system have a geometrical interpretation. 
Let us define the extended order vectors a® = (l,a) S and extended 

exponent vectors a® = (7^,0:^) G Let us consider the equality pLj = pj/. 

This represents the equation of a n dimensional hyperplane of ]R"+^, orthogonal 
to the vector ccj® — ccj'®: 

(a®,a/) = (a®,aj/®), (7) 

where (,) is the dot product in We will see in the next section that 

the minimality condition on the exponents pLj implies that the normal vectors 
aj® — aj'® are edges of the so-called Newton polytope [aim. 

For each equation i, let us define 

Mi{a) = argmin(/rj(a), S'jj > 0 ) = argmin(/rj(a), < 0 ), (8) 

3 3 

in other words Mi denote the set of monomials having the same minimal expo¬ 
nent /ii. 

We call tropically truncated system, the system obtained by keeping only the 
dominating monomials in Q, as follows: 

^ ( 9 ) 

jGMi(a) 

The tropical truncated system is uniquely determined by the index sets 
Mi{a), therefore, by the tropical equilibration a. Reciprocally, two tropical 
equilibrations can have the same index sets Mi{a) and truncated systems. 
We say that two tropical equilibrations ai, 02 are equivalent iff Mi{ai) = 
Mi{a2), for all i. Equivalence classes of tropical equilibrations are called branches. 
For each branch there exists a unique convex polytope, cf. Q. The union of 
branches are subsets of tropical prevariety. It is a subset because we are in¬ 
terested in tropical equilibration of at least two monomials of different signs 
as expressed in Q. This sign condition is essential as we are interested to 
approximate positive real solutions of the polynomial system. 
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Minimal Branches A branch B with an index set Mi is minimal if M' C Mi 
for all i where M/ is the index set of a branch B' implies B' = B or B' = 
emptyset. In the terminology of convex polytopes, this means a branch B 
with a convex polytope Pi is minimal if Pi C P[ for all i where P/ is the convex 
polytope for branch B' implies B' = B or B' is empty. For each index i, relation 
0 defines a hyperplane, the tropical equilibration branches are on intersections 
of k such hyperplanes where k is number of polynomial equations representing 
the right hand side of Q. Minimal branches are maximal (w.r.t. inclusion) 
polytopes in the tropical prevariety. 

Connected Components of Minimal Branches Two minimal branches 
represented by index sets Mi and Mj are connected if there exists a branch 
with index set Mk such that Mi C Mk and Mj C Mk- In the terminology of 
convex polytopes, this amounts to checking the intersection between two convex 
polytopes Pi and Pj (corresponding to minimal branches Mi and Mj) if whether 
Pi n Pj is non void for all i ^ j. 


4 Algorithm 

In this section, we describe an algorithm to compute the tropical equilibrations, 
test the equilibrations for the equivalence classes (i.e., branches) and compute 
the minimal branches (cf. section]^. 

4.1 Newton Polytope and Edge Filtering 

Given the input polynomial in the form of 0. for each equation and species i, 
we define a Newton polytope A/), that is the convex hull of the set of points a® 
such that Sij ^ 0 and also including together with all the points the half-line 
emanating from these points in the positive e direction. This is the Newton 
polytope of the polynomial in right hand side of 0 , with the scaling parameter 
e considered as a new variable. 

As explained in section the tropical equilibrations correspond to vectors 
a® = (l,a) S satisfying the optimality condition as per ([^. This con¬ 

dition is satisfied automatically on hyperplanes orthogonal to edges of Newton 
polytope connecting vertices a®,, ctj,, satisfying the opposite sign condition. 
Therefore, a subset of edges from the Newton polytope is selected based on the 
filtering criteria which tells that the vertices belonging to an edge should be 
from opposite sign monomials as explained in ( |I0[ ). 

E{P) = {{vi,V2} C (J) I conv(?;i,?;2) G Fi{P) 

Asign(ui) X sign(r; 2 ) = -1}, (10) 

where Vi is the vertex and V is the vertex set of the Newton polytope, conv(r;i, V 2 ) 
is the convex hull of vertices vi,V 2 and Fi{P) is the set of 1-dimensional face 
(edges) of the Newton polytope, sign(ui) represents the sign of the monomial 
which corresponds to vertex Vi. Figureshows an example of Newton polytope 
construction for a single equation. Further definitions about properties of a 
polytope and Newton polytope can be found in 
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Figure 1: An example of a Newton polytope for the polynomial —xf+XiX 2 —xl+XiX 2 - 
In this example, all the monomial coefficients have order zero in e and we want to solve 
the tropical problem min(3ai + a 2 , ai + 202 ) = min(6ai, 3ai). The Newton polytope 
vertices (6, 0), (3,0), (1, 2) are connected by lines. The point (3,1) is not a vertex as 
it lies in the interior of the polytope. This stems to having min(3ai + 02 , ai + 202 ) = 
oi + 202 for all tropical solutions, which reduces the number of cases to be tested. 
The thick edges satisfy the sign condition, whereas the dashed edge does not satisfy 
this condition. For this example, the solutions of the tropical problem are in infinite 
number and are carried by the two half-lines oi = 02 > 0 and 5ai = 2 a 2 < 0, 
orthogonal to the thick edges of the Newton polygon 


4.2 Computing Branches of Tropical Equilibrations 

Using the Newton polytope formulation, one can then solve the tropical equili¬ 
bration problem in ^ using the edges of Newton polytope (as in Q). A feasible 
solution is a vector (oi,... ,a„) satisfying all the equations of system (§ and 
lies in the intersection of hyperplanes (or convex subsets of these hyperplanes) 
orthogonal to edges of Newton polytopes obeying the sign conditions. Of course, 
not all sequences of edges lead to non-void intersections and, thus, feasible solu¬ 
tions. This can be tested by the following linear programming problem resulting 
from ([^: 

7j(f) -I- {a,aj{i)) = ij(i) + (a,Q;'.(i)) < 7" -h 

for allj" i = l,...,n (11) 

where j(i),/(i) define the chosen edge of the fth Newton polytope. The set 
of indices j" can be restricted to vertices of the Newton polytope, because the 
inequalities are automatically fulfilled for monomials that are internal to the 
Newton polytope. From ( |IT| ) , the sequence of edges leading to a feasible solution 
is actually a set of linear inequalities and hence constitutes a feasible solution 
system (convex polytope), which was computed using Algorithm]^ Such feasible 
solution systems are actually convex polytopes as defined in §. The equivalence 
classes among the feasible solution systems constitute the equivalence classes 
of tropical equilibrations called branches. This was done using the method 
equaLpolyhedra implemented in the software package polymake [5]. For instance, 
in the example of the preceding section, the choice of the edge connecting vertices 
(1, 2) and (3, 0) leads to the following linear programming problem: 

oi + 202 = 3oi < 6ai, 
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whose solution is a half-line orthogonal to the edge of the Newton polygon. The 
pseudo-code is presented in Algorithm It is clear from the above that the 
possible choices are exponential. In order to improve the running time of the 
algorithm, the pruning strategy evaluates ( [Tl] ) in several steps (cf. Algorithm 
and Fig. [^. It starts with an arbitrary pair of edges and proceeds to add the 
next edge only when the inequalities (11) restricted to these two pair of edges 
are satisfied. The pruning method is a heuristic to filter out the infeasible set 
of edge combinations. A similar approach was undertaken in 



Figure 2 : Pruning strategy. The possible combinations of edges are represented in a 
tree representation where aj represents ith edge from j'th Newton polytope. An edge 
set nci is the set of edges for ith Newton polytope. The algorithm starts by testing for 
feasible solution for first pair of edge sets. If a feasible solution is found, the algorithm 
proceeds further to other edge sets or it backtracks. In the figure, en and 621 are 
selected from edge sets nei, ne2 and are checked for a feasible solution satisfying O- 
If such a solution exists, it moves to esifrom the next edge set and again checks for 
feasible solution, if not then it backtracks to 621 and then to 632 which results in a 
feasible solution. Therefore, the sub-tree with root node 631 is discarded from future 
searches and this improves running time. Likewise the branch en and 622 is explored. 
This approach is similar to the branch and bound algorithm technique. The dashed 
arrows show the flow of the program 


4.3 Computation of Minimal Branches 

The minimal branches are explained in section Computation was performed 
using included-polyhedra method in polymake. The minimal branches as well as 
branches contained in minimal branches are represented in Fig. [^as a directed 
graph with layers where top most layer are minimal branches. 

4.3.1 Sample Point for Minimal Branches 

The polytopes corresponding to minimal solution branches obtained from the 
previous steps are represented by their facets and affine hull properties which are 
basically the set of inequalities and equalities. From such a set of inequalities, 
a sample point (ai,...,a„) is computed using Satisfiability Modulo Theories 
(SMT) solver called Microsoft Z3 software [5] in python programming environ¬ 
ment. With Microsoft Z3, one can generate the sample point belonging exclu- 
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Algorithm 1 : SolveOrders: Steps of tropical equilibration algorithm 


Input: List of edge sets nei,ne 2 , ...,ne„ (cf. Fig. [^, and the 
corresponding vertices of Newton polytope 
Output: Set of feasible solution systems (convex polytopes) 
corresponding to orders of the variables Oi, 02 ,a„ 
equilibration solution set) 


1 begin 

2 solutionset ={}; integer fc=l; equation = {} 
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SolveOrders(equation, k, edge-sets, vertices) 


4 

5 

6 

7 

8 


if k > n then 

1^ return 


for I = 1 to number of entries in nekedge-set do 
equation(k)* = vertices in /*^row 

inequalities* = all other vertices in neito ne^edge-sets 
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10 

11 

12 


if LinearSolve(equation,inequalities)is feasible then 
if k = n then 

^ add the feasible solution system to solutionset 
SolveOrders(equation, fc -I- 1, nei, ..^nek, vertices) 


(tropical 


13 *The equations and inequalities are initialised as per 0 


sively to a minimal branch (and not at the intersections of minimal branches) 
by corresponding Boolean conjugations as shown in the following manner 

{T,€B,AT,iN„yiGl} (12) 

where Ti is a tropical equilibration solution corresponding to Bi where B is 
set of polytopes corresponding to minimal solutions and N is the set of rest 
solution branches along with minimal solution branches B \ Bi . I is an index 
set denoting the elements of B. 

The sample point thus obtained is a feasible solution to 0. For our pur¬ 
pose, the benefit of using Z3 over any existing linear programming software is 
that it distinguishes strict and non-strict inequality conditions, which allows us 
to generate the sample point belonging exclusively to a minimal branch. 

5 Results 

To compute the tropical equilibrations and to demonstrate the running time 
of our algorithm, 33 models from the r25 version of Biomodels database m 
having polynomial vector field were parsed with PoCaB m- 

5.1 Summary 

A summary of the analysis is presented in Table The analysis is performed 
to compute all possible combinations of Newton polytope edges leading to trop¬ 
ical solutions within a maximal running time of 10000 seconds of CPU time. 







In practice, we restrict this search space using the tree pruning strategy as ex¬ 
plained in Fig. The analysis was repeated with four different choices for 
e values. In our framework, the number of variables may not be equal to the 
number of equations as the conservation laws (that are sums of variables whose 
total concentration is invariant) are treated as extra linear equations in our 
framework. 

5.2 Running Time 

A semilog time-plot is presented in Fig. Qa) which plots the log of running time 
in milliseconds versus the number of equations in the model. In Fig. |^b), the 
pruning ratio, i.e., the efficacy of tree pruning for e value of 1/5 is plotted. The 
pruning ratio is the ratio between the number of times the linear programming 
is invoked with every tree pruning step (cf. Fig. and the possible number 
of combinations of Newton polytope edges possible without tree pruning. This 
ratio is, thus, a measure of efficiency achieved due to pruning. 



BIOMD0000000035 BIOMD0000000072 

Figure 3: A directed graph in layered form showing the inclusion relations among the 
different solution branches (for e = 1/11) for two models namely BIOMDOOOOOOOO- 
35,72. Vertices in the graph comprise of polytopes corresponding to solution branches 
and an directed edge between i and j means j is included in i. The topmost layer 
contain the minimal solution branches, thereafter the bottom layers are ’’included” 
solution branches. The layers of the included solution branches are based on the 
dimension of the corresponding pol 3 rtopes (arranged in descending order). Therefore, 
included solutions in one layer are of same dimensional polytope 


5.3 Minimal Branches and Role of e 

A semilog plot for minimal solution branches is presented in Fig. j^a) and a 
semilog plot in Fig. [^b) showing the ratio of minimal solution branches to the 
number of feasible solution systems (obtained from Algorithm [^. It shows that 
a large proportion of feasible solution systems are either redundant or included 
in other feasible systems (i.e., inclusion relations). 
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Table 1: Summary of analysis on Biomodels database. Tropical solutions here mean 
existence of at least one feasible solution from all possible combination of vertices of the 
Newton polytope. Timed-out means all solutions could not be computed within 10000 
secs of computation time. No tropical solution implies that no possible combination of 
edges could be found resulting in a feasible solution. Model BIOMD0000000289 has so¬ 
lution at £ values 1/5,1/7 and 1/9 but no solutions at 1/23. Model BIOMD0000000108 
has no solutions at all values of £ 


e 

value 

Total 

models 

con¬ 

sidered 

Models 

without 

tropical 

solutions 

Models 

with 

tropical 

solutions 

Average 
running 
time (in 
secs) 

Average 
number of 

minimal 

branches 

1/5 

33 

1 

32 

299.38 

3.24 

1/7 

33 

1 

32 

244.12 

3 

1/9 

33 

1 

32 

309.73 

3.75 

1/23 

33 

2 

31 

3179.32 

3.84 
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Figure 4: (a) Plot of running time against number of equations in the model, (b) 

Pruning ratio for e value of 1/5 against number of equations in the model 


In order to investigate the effect of different e values on the number of 
minimal solutions, a boxplot is presented in Fig. [^a) for different choices of e 
values. In Fig. ib), the boxplot shows the ratio of minimal solution branches 
to the number of feasible solution systems (obtained from Algorithm for 
different choices of e values 
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10 20 30 40 10 20 30 40 

Number of equations Number of equations 


Figure 5: (a) Minimal branches against number of equations in the model, (b) Ratio 

of minimal branches to the number of feasible solution systems, i.e., number of feasible 
edge combinations against number of equations in the model 


5.4 Slow-Fast Variables 

From the tropical equilibrations, we computed fj,j — ai{ci. @) which allow us 
to order the variables of the model, from the fastest (smallest — a^) to the 
slowest (largest — ai). This is a measure to separate the variables into slow 
and fast which is an important step in constructing the invariant manifold for 
model reduction [Ml m [la E]. The heatmap in Fig. [^and Fig. plots this 
value for minimal solution branch and all solution branches, respectively, for 
few selected models. For some models, there appears to be a natural clustering 
(as seen from the dendogram from the hierarchical clustering) which requires 
further investigation. 

5.5 Connected Components 

As described in section and defined in section tropical solutions can be 
roughly associated with metastable states and these branches (which are convex 
polytopes) form a polyhedral complex. The number of connected components of 
this complex indicates the possible number of transitions between the minimal 
solution branches. Figure depicts such a graph of connected components of 
minimal solutions branches. 
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Figure 6: Boxplots showing (a) Distribution of Minimal branches, (b) Ratio of 
minimal branches to the number of feasible solution systems. Both distributions are 
at different e values: 1/5,1/7,1/9,1/23 



BIOMD0000000035 


BIOMD0000000072 


Figure 7: Heatmaps showing the rescaled orders (at e = 1/11) for four models 
namely BIOMDOOOOOOOO-35,72 with hierarchical clustering for variables (horizontal 
axis) and tropical solutions (vertical axis). 


6 Discussions 

We present an algorithm to compute the tropical equilibrations and organise 
them into branches and minimal branches. The directed graphs showing the 
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BIOMD0000000035 BIOMD0000000072 

Figure 8: Heatmaps showing the rescaled orders (at e = 1/11) for four models 
namely BIOMDOOOOOOOO-35,72 with hierarchical clustering for variables (horizontal 
axis) and tropical solutions (vertical axis). The heatmaps include the minimal solution 
branches (with prehx ” M”) and other solution branches excluding minimal ones (with 
prehx ”S”) 



BIOMD0000000035 BIOMD0000000072 

Figure 9: Graph of connected components (at e = 1/11) for four models namely 
BIOMDOOOOOOOO-35,72. All of them have one connected component. The vertices are 
minimal solution branch and there exists an edge if the intersection between the two 
vertices is non-void 


inclusion relations between branches of tropical equilibration solutions reveal a 
rich structure. In addition, the connectivity of minimal branches are computed 
which provides an estimate for the possible dynamical transitions between them. 
One of the applications of tropical equilibration in systems biology is model 
reduction as demonstrated in [IllIT]. The concentration orders depicted in the 
heatmaps demonstrate the applicability of the algorithm in this direction. Thus, 
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the overall solution structure provides insights into the dynamics of the system. 
More precisely, for the biochemical models in the Biomodels database, a large 
number of feasible solution systems were obtained but the number of minimal 
branches is considerably less. For example, the number of minimal branches at 
e = 1/5 ranged from 1 to 17, whereas the total number of solutions ranged from 
1 to 9847. 

As the dominant terms in the polynomial system are the same for all the 
tropical solution on branches, it could be that branches correspond to invariant 
manifolds. In the same spirit, minimal branches could correspond to minimal 
invariant manifolds. This idea will be pursued in future work. Lastly, we have 
shown an application of tropical geometry to invariant manifolds defined by 
polynomial systems but the direct application of tropical geometry to differential 
equation systems mm is also known. 
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